module objects.triangle;

import utils.vector3;
import utils.vector2;
import scene.ray;
import scene.intersection;
import std.math;

static const float DIV_EPSILON = (1e-10f); // Avoid numerically instable divisions

int g_IntersectCount = 0;

version = MedianPoint;

struct Triangle {
    Vector3f[3] v;

    /**
     * Intersection of a ray with a triangle
     * Algorithm from paper
     * Fast Minimum Storage Ray Triangle Intersection
     * Tomas Moeller, Ben Trumbore
     */
    public bool intersectRay(inout Ray isRay, inout Intersection intersect) {
        static Vector3f e1, e2, tvec, pvec, qvec;
        static float det, invDet;

        debug g_IntersectCount++;

        //sub3f(e1,v[1],v[0])
        e1.x1 = v[1].x1 - v[0].x1;
        e1.x2 = v[1].x2 - v[0].x2;
        e1.x3 = v[1].x3 - v[0].x3;

        //sub3f(e2,v[2],v[0])
        e2.x1 = v[2].x1 - v[0].x1;
        e2.x2 = v[2].x2 - v[0].x2;
        e2.x3 = v[2].x3 - v[0].x3;

        //cross3f(pvec,isRay.dir,e2)
        pvec.x1 = isRay.dir.x2*e2.x3 - isRay.dir.x3*e2.x2;
        pvec.x2 = isRay.dir.x3*e2.x1 - isRay.dir.x1*e2.x3;
        pvec.x3 = isRay.dir.x1*e2.x2 - isRay.dir.x2*e2.x1;

        //det = e1 * pvec
        det = e1.x1*pvec.x1 + e1.x2*pvec.x2 + e1.x3*pvec.x3;

        version(NoCull) {
            /* the non-culling branch */

            if (det > -DIV_EPSILON && det < DIV_EPSILON)
                return false;

            invDet = 1.0f / det;

            //sub3f(tvec,isRay.start,v[0])
            tvec.x1 = isRay.start.x1 - v[0].x1;
            tvec.x2 = isRay.start.x2 - v[0].x2;
            tvec.x3 = isRay.start.x3 - v[0].x3;

            //intersect.u = (tvec * pvec) * invDet
            intersect.u = tvec.x1*pvec.x1 + tvec.x2*pvec.x2 + tvec.x3*pvec.x3;
            intersect.u *= invDet;
            if (intersect.u < 0.0f || intersect.u > 1.0f)
                return false;

            //cross3f(qvec,tvec,e1)
            qvec.x1 = tvec.x2*e1.x3 - tvec.x3*e1.x2;
            qvec.x2 = tvec.x3*e1.x1 - tvec.x1*e1.x3;
            qvec.x3 = tvec.x1*e1.x2 - tvec.x2*e1.x1;

            //intersect.v = (isRay.dir * qvec) * invDet
            intersect.v = isRay.dir.x1*qvec.x1 + isRay.dir.x2*qvec.x2 + isRay.dir.x3*qvec.x3;
            intersect.v *= invDet;
            if (intersect.v < 0.0f || intersect.u + intersect.v > 1.0f)
                return false;

            //intersect.t = (e2 * qvec) * invDet
            intersect.t = e2.x1*qvec.x1 + e2.x2*qvec.x2 + e2.x3*qvec.x3;
            intersect.t *= invDet;
        } else {
            /* culling branch */

            if (det < DIV_EPSILON)
                return false;

            //sub3f(tvec,isRay.start,v[0])
            tvec.x1 = isRay.start.x1 - v[0].x1;
            tvec.x2 = isRay.start.x2 - v[0].x2;
            tvec.x3 = isRay.start.x3 - v[0].x3;

            //intersect.u = tvec * pvec
            intersect.u = tvec.x1*pvec.x1 + tvec.x2*pvec.x2 + tvec.x3*pvec.x3;
            if (intersect.u < 0.0f || intersect.u > det)
                return false;

            //cross3f(qvec,tvec,e1)
            qvec.x1 = tvec.x2*e1.x3 - tvec.x3*e1.x2;
            qvec.x2 = tvec.x3*e1.x1 - tvec.x1*e1.x3;
            qvec.x3 = tvec.x1*e1.x2 - tvec.x2*e1.x1;

            //intersect.v = isRay.dir * qvec
            intersect.v = isRay.dir.x1*qvec.x1 + isRay.dir.x2*qvec.x2 + isRay.dir.x3*qvec.x3;
            if (intersect.v < 0.0f || intersect.u + intersect.v > det)
                return false;

            //intersect.t = e2 * qvec
            intersect.t = e2.x1*qvec.x1 + e2.x2*qvec.x2 + e2.x3*qvec.x3;
            invDet = 1.0f / det;
            intersect.t *= invDet;
            intersect.u *= invDet;
            intersect.v *= invDet;
        }
        return true;
    }

    public static Triangle opCall(Vector3f[3] v) {
        Triangle ret;
        ret.v[] = v[];
        return ret;
    }

    ///get triangle normal
    public Vector3f normal() {
        return ((v[1] - v[0])|(v[2] - v[0])).normal();
    }

    public char[] toString() {
        return v[0].toString() ~ ", " ~ v[1].toString() ~ " " ~ v[2].toString();
    }

    public Vector3f medianPoint() {
        return (v[0]+v[1]+v[2])/3;
    }

    float getMin(uint axis) {
        float m = v[0][axis];
        float a = v[1][axis];
        float b = v[2][axis];
        if (a < m)
            m = a;
        if (b < m)
            m = b;
        return m;
    }

    float getMax(uint axis) {
        float m = v[0][axis];
        float a = v[1][axis];
        float b = v[2][axis];
        if (a > m)
            m = a;
        if (b > m)
            m = b;
        return m;
    }

    /** Check if triangle belongs to "left" or "right" side
     * leftPlane and rightPlane are adapted accordingly and will contain the
     * actual splitplanes after all tris have been processed
     * Returns: true if triangle belongs to left group
     */
    bool isLeft(float split, int axis, inout float leftPlane, inout float rightPlane) {
        auto min = getMin(axis);
        auto max = getMax(axis);

        assert(min <= max);

        if (max <= split) {
            //object on left side
            if (max > leftPlane)
                leftPlane = max;
            return true;
        } else if (min >= split) {
            //object on the right side
            if (min < rightPlane)
                rightPlane = min;
            return false;
        } else {
            //object crosses split plane, put it to where there is more
            version(MedianPoint) {
                float m = (v[0][axis]+v[1][axis]+v[2][axis])/3;
                bool left = m < split;
            } else {
                bool left = split-min > max-split;
            }
            if (left) {
                if (max > leftPlane)
                    leftPlane = max;
                return true;
            } else {
                if (min < rightPlane)
                    rightPlane = min;
                return false;
            }
        }
        assert(false); //never reached
    }
}

const uint modulo[] = [0,1,2,0,1];

struct TriAccel {
    // first 16 byte half cache line
    // plane:
    float n_u; //!< == normal.u / normal.k
    float n_v; //!< == normal.v / normal.k
    float n_d; //!< constant of plane equation
    int k; // projection dimension
    // second 16 byte half cache line
    // line equation for line ac
    float b_nu;
    float b_nv;
    float b_d;
    //int pad1; // pad to next cache line
    // third 16 byte half cache line
    // line equation for line ab
    float c_nu;
    float c_nv;
    float c_d;
    //int pad2; // pad to 48 bytes for cache alignment purposes

    public void calculate(inout Triangle tri) {
        Vector3f b = tri.v[2] - tri.v[0];
        Vector3f c = tri.v[1] - tri.v[0];
        Vector3f N = c | b;

        // determine projection dimensiondimensions
        if (abs(N.x) > abs(N.y))
            if (abs(N.x) > abs(N.z)) k = 0; /* X */ else k=2; /* Z */
        else
            if (abs(N.y) > abs(N.z)) k = 1; /* Y */ else k=2; /* Z */
        int u = (k+1) % 3; int v = (k+2) % 3;

        Vector3f Ns = N/N[k];

        n_u = Ns[u];
        n_v = Ns[v];
        n_d = tri.v[0] * Ns;

        b_nu = -b[v]/(b[u]*c[v] - b[v]*c[u]);
        b_nv = b[u]/(b[u]*c[v] - b[v]*c[u]);
        b_d = (b[v]*tri.v[0][u] - b[u]*tri.v[0][v])/(b[u]*c[v] - b[v]*c[u]);

        c_nu = c[v]/(b[u]*c[v] - b[v]*c[u]);
        c_nv = -c[u]/(b[u]*c[v] - b[v]*c[u]);
        c_d = (-c[v]*tri.v[0][u] + c[u]*tri.v[0][v])/(b[u]*c[v] - b[v]*c[u]);
    }

    bool intersectRay(inout Ray ray, inout Intersection hit) {
        debug g_IntersectCount++;

        // start high-latency division as early as possible
        int u = modulo[k+1]; int v = modulo[k+2];
        float nd = 1./(ray.dir[k]
            + n_u * ray.dir[u] + n_v * ray.dir[v]);
        float f = (n_d - ray.start[k]
            - n_u * ray.start[u] - n_v * ray.start[v]) * nd;

        // check for valid distance.
        if (f <= 0.0f || f >= ray.tMax)
            return false;

        // compute hitpoint positions on uv plane
        float hu = (ray.start[u] + f * ray.dir[u]);
        float hv = (ray.start[v] + f * ray.dir[v]);

        // check first barycentric coordinate
        float lambda = (hu * b_nu + hv * b_nv + b_d);
        if (lambda < 0.0f)
            return false;

        // check second barycentric coordinate
        float mue = (hu * c_nu + hv * c_nv + c_d);
        if (mue < 0.0f)
            return false;

        // check third barycentric coordinate
        if (lambda+mue > 1.0f)
            return false;

        // have a valid hitpoint here. store it.
        hit.t = f;
        hit.u = lambda;
        hit.v = mue;
        return true;
    }

};

struct TriangleNorms {
    Vector3f[3] n;

    public static TriangleNorms opCall(Vector3f n1, Vector3f n2, Vector3f n3) {
        TriangleNorms ret;
        ret.n[0] = n1;
        ret.n[1] = n2;
        ret.n[2] = n3;
        return ret;
    }

    public static TriangleNorms opCall(Vector3f[3] n) {
        TriangleNorms ret;
        ret.n[] = n;
        return ret;
    }

    public Vector3f getNormal(float u, float v) {
        return (n[0]*(1.0f-u-v)+n[1]*u+n[2]*v).normal();
    }

    public char[] toString() {
        return n[0].toString() ~ ", " ~ n[1].toString() ~ " " ~ n[2].toString();
    }
}

struct TriangleTex {
    Vector2f[3] texCoords;

    public static TriangleTex opCall(Vector2f[3] texCoords) {
        TriangleTex ret;
        ret.texCoords[] = texCoords[];
        return ret;
    }
}
